# x: the vector
# n: the number of samples
# centered: if FALSE, then average current sample and previous (n-1) samples
#           if TRUE, then average symmetrically in past and future. (If n is even, use one more sample from future.)
movingAverage <- function(x, n=1, centered=FALSE) {

    if (centered) {
        before <- floor  ((n-1)/2)
        after  <- ceiling((n-1)/2)
    } else {
        before <- n-1
        after  <- 0
    }

    # Track the sum and count of number of non-NA items
    s     <- rep(0, length(x))
    count <- rep(0, length(x))

    # Add the centered data 
    new <- x
    # Add to count list wherever there isn't a 
    count <- count + !is.na(new)
    # Now replace NA_s with 0_s and add to total
    new[is.na(new)] <- 0
    s <- s + new

    # Add the data from before
    i <- 1
    while (i <= before) {
        # This is the vector with offset values to add
        new   <- c(rep(NA, i), x[1:(length(x)-i)])

        count <- count + !is.na(new)
        new[is.na(new)] <- 0
        s <- s + new

        i <- i+1
    }

    # Add the data from after
    i <- 1
    while (i <= after) {
        # This is the vector with offset values to add
        new   <- c(x[(i+1):length(x)], rep(NA, i))

        count <- count + !is.na(new)
        new[is.na(new)] <- 0
        s <- s + new

        i <- i+1
    }

    # return sum divided by count
    s/count
}

deviation <- function(v1, v2){

    i <- 1
    count <- min(length(v1), length(v2))
    deviation <- 0
    for(i in (1:count)){
        deviation <- abs(v1[i]^2 - v2[i]^2)
    }
    deviation

}

#set current directory
setwd("C:/Project/Development/nvrocket/R")

#load data
first = read.csv("MF+T1+SDB.csv")
second = read.csv("SH+T1+SDB.csv")

attach(first)
attach(second)

# set directory
plot(first$Ping_time, first$Depth)
plot(second$Ping_time, second$Depth)

# find maximal covaration

#find longest vector
firstLen = length(first$Depth)
secondLen = length(second$Depth)
minLen = min(firstLen, secondLen)
maxLen = max(firstLen, secondLen)
minVector = movingAverage(first$Depth, 21, TRUE)
maxVector = movingAverage(second$Depth, 21, TRUE)
if(firstLen > secondLen){
    minVector = movingAverage(second$Depth, 21, TRUE)
    maxVector = movingAverage(first$Depth, 21, TRUE)
}

minSelectionLen = 300
maxCor = 0
dev = 999999
curDev = -1
positionMaxStart = -1
positionMaxEnd = -1
positionMinStart = -1
positionMinEnd = -1
#calculate corelation from the start
for(i in (minSelectionLen:minLen)){
    startIndex = minLen - i + 1
    curCor = abs( cor(minVector[startIndex : minLen], maxVector[1 : i], method = "pearson") )
    if(curCor > maxCor){
        maxCor = curCor
#       positionMaxStart = 1
#       positionMaxEnd = i
#       positionMinStart = startIndex
#       positionMinEnd = minLen
    }
    curDev <- deviation(minVector[startIndex : minLen], maxVector[1 : i])
    if(curDev < dev){
         positionMaxStart = 1
         positionMaxEnd = i
         positionMinStart = startIndex
         positionMinEnd = minLen    

         dev = curDev
    }
}
#calculate in the middle
for( i in (1: (maxLen - minLen) ) ){
    startMinIndex = 1;
    endMinIndex = minLen;
    startMaxIndex = i;
    endMaxIndex = startMaxIndex + minLen - 1;
    curCor = abs( cor(minVector[startMinIndex : endMinIndex], maxVector[startMaxIndex : endMaxIndex], method = "pearson") )
    if(curCor > maxCor){
        maxCor = curCor
#        positionMaxStart = startMinIndex
#        positionMaxEnd = endMinIndex
#        positionMinStart = startMaxIndex
#        positionMinEnd = endMaxIndex       
    }
    curDev <- deviation(minVector[startMinIndex : endMinIndex], maxVector[startMaxIndex : endMaxIndex])
    if(curDev < dev){
         positionMaxStart = 1
         positionMaxEnd = i
         positionMinStart = startIndex
         positionMinEnd = minLen    

         dev = curDev
    }    
}
#calculate corelation to the end
for(i in (minSelectionLen:minLen)){
    startMaxIndex = maxLen - i + 1
    curCor = abs( cor(minVector[1 : i], maxVector[startMaxIndex : maxLen], method = "pearson") )
    if(curCor > maxCor){
        maxCor = curCor
        positionMaxStart = startMaxIndex
        positionMaxEnd = maxLen
        positionMinStart = 1
        positionMinEnd = i
    }
    curDev <- deviation(minVector[1 : i], maxVector[startMaxIndex : maxLen])
    if(curDev < dev){
         positionMaxStart = 1
         positionMaxEnd = i
         positionMinStart = startIndex
         positionMinEnd = minLen    

         dev = curDev
    }  
    
}

#hardcoded
plot(first$Ping_time[positionMaxStart:positionMaxEnd], first$Depth[positionMaxStart:positionMaxEnd])
plot(second$Ping_time[positionMinStart:positionMinEnd], second$Depth[positionMinStart:positionMinEnd])


# clear resources
detach(first)
detach(second)
